function [imgA] = gencamimg_plane(intrinsics, imgW, imRefW, CWC, rOCOWC, noiseStdDev, fillValues, fileNamePrefix, fileExt, fileName)
%GENCAMIMG_PLANE 生成世界坐标系下的平面图像在相机中的成像图像
%
% Tips
% # 当使用imwarp函数应用本函数输出的位移场D时，为获得正确的畸变，需要取负值，即imwarp(I, -D)
%
% Input Arguments
% # intrinsics: cameraIntrinsics对象，相机内参
% # imgW: 矩阵，世界坐标系W（与平面固连，一般定义XW与YW轴分别沿imgW图像的横向及纵向）下的平面灰度图像
% # imRefW: imref2d对象，描述平面图像内在坐标（intrinsic coordinates）与世界坐标系坐标的关系
% # CWC: 3*3*p数组，n为需要生成的图像张数，每张图像对应的世界坐标系W相对于相机坐标系C的姿态
% # rOCOWC: 3*p矩阵，每张图像对应的世界坐标系W原点相对于相机坐标系C原点的位矢
% # noiseStdDev: 标量，图像噪声（详见distortimage函数），空或者NaN表示不加噪声，单位像素，默认为空
% # fillValues: 标量，畸变后图像边界未赋值部分的填充值，默认为255（白色）
% # fileNamePrefix: 字符串，图片文件前缀（含路径），图片名称为“前缀数字编号.扩展名”，默认为空
% # fileExt: 字符串，图片文件扩展名，当fileNamePrefix及fileExt全为空时不保存图片文件，默认为png
% # fileName: p*1字符元包数组，每张图像对应的图片文件名，默认为空
%
% Output Arguments
% # imgA: 与imgW类型相同的m*n*p矩阵，其中m、n分别为intrinsics.ImageSize的第2、1个元素，平面在相机中成像的图像，考虑单应变换、畸变及噪声

if nargin < 6
    noiseStdDev = [];
end
if nargin < 7
    fillValues = 255;
end
if nargin < 8
    fileNamePrefix = '';
end
if nargin < 9
    fileExt = 'png';
end
if nargin < 10
    fileName = [];
end

nImg = size(CWC, 3);
assert(size(rOCOWC, 2) == nImg);

xResolution = intrinsics.ImageSize(1);
yResolution = intrinsics.ImageSize(2);

if nargout > 0
    imgA = NaN(size([yResolution, xResolution nImg]));
else
    imgA = [];
end

outputIsNeeded = ~isempty(imgA);
parfor i=1:nImg
    % project onto camera imaging plane
    H = intrinsics.IntrinsicMatrix'*[CWC(:, 1, i) CWC(:, 2, i) rOCOWC(:, i)]; % homography matrix
    [thisImgA_full, imgRefA_full] = imwarp(imgW, imRefW, projective2d(H'), 'FillValues', fillValues);
    
    % locate camera field of view
    [xIntr, yIntr] = worldToIntrinsic(imgRefA_full, [0.5 xResolution+0.5], [0.5 yResolution+0.5]);
    xIntr = round(xIntr);
    yIntr = round(yIntr);
    thisImgA_FOV = fillValues * ones(diff(yIntr)+1, diff(xIntr)+1, class(thisImgA_full)); % in the field of view of the camera
    xInd = [max(xIntr(1), 1) min(xIntr(2), imgRefA_full.ImageSize(2))];
    yInd = [max(yIntr(1), 1) min(yIntr(2), imgRefA_full.ImageSize(1))];
    thisImgA_FOV((yInd(1):yInd(2))+(1-yIntr(1)), (xInd(1):xInd(2))+(1-xIntr(1))) = thisImgA_full(yInd(1):yInd(2), xInd(1):xInd(2));
    thisImgA_full = []; % clear to save memory
    
    if isNaN(noiseStdDev) || isempty(noiseStdDev)
        thisImgA = imresize(thisImgA_FOV, [yResolution, xResolution], 'bilinear');
    else
        % add distortion, noise, and resize to fit actual camera image size
        scale = fliplr(size(thisImgA_FOV)) ./ intrinsics.ImageSize;
        intrinsics_FOV = cameraIntrinsics(intrinsics.FocalLength .* scale, intrinsics.PrincipalPoint .* scale, fliplr(size(thisImgA_FOV)), ...
            'Skew', intrinsics.Skew, 'RadialDistortion', intrinsics.RadialDistortion, 'TangentialDistortion', intrinsics.TangentialDistortion);
        distimage = distortimage(intrinsics_FOV, thisImgA_FOV, mean(scale)*noiseStdDev, fillValues);
        thisImgA = imresize(distimage, [yResolution, xResolution], 'bilinear');
    end
    
    % ouput
    if outputIsNeeded
        imgA(:, :, i) = thisImgA;
    end
    if ~isempty(fileName)
        imwrite(thisImgA, fileName{i, 1});
    elseif ~isempty(fileNamePrefix) || ~isempty(fileExt)
        imwrite(thisImgA, [fileNamePrefix num2str(i), '.', fileExt]);
    end
end
end